#############################################################
#################       Linear Polynomial        #############
##############   Robustness on Education 
############  December 15, 2022 

rm(list=ls())
library(foreign)
library(plyr)
library(readstata13)
library(multiwayvcov)
library(sandwich)
library(lmtest)
library(stargazer)


data=read.csv("~/Dropbox/Personal Research 2017/replications/karn_nov16.csv")
names(data)
summary(data$Latitude)
summary(data$Longitude)
summary(data$border1)
summary(data$border2)
summary(data$Slope)


#######################
##### Distances #####

# 20 km
#Distance to Mysore-Bombay Border
rd10.mb=data[which(data$NEAR_DIST_border1<10000),] #20 km

table(rd10.mb$border1)

#Distance to Hyderabad-Bombay Border
rd10.hb=data[which(data$NEAR_DIST_border2<10000),] #20 km

table(rd10.hb$border2)


#baseline bandwidth
#Mysore
#outcome-high schools
educ.mys10=lm(high_binary~border1+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd10.mb) #OLS estimation
summary(educ.mys10)
educ.mys10.cl=cluster.vcov(educ.mys10, rd10.mb$dist_name)
educ.mys10.se=sqrt(diag(educ.mys10.cl)) #cluster standard errors


#Hyderabad
#outcome-high schools
educ.hyd10=lm(high_binary~border2+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd10.hb) #ols estimation 
summary(educ.hyd10)
educ.hyd10.cl=cluster.vcov(educ.hyd10, rd10.hb$dist_name)
educ.hyd10.se=sqrt(diag(educ.hyd10.cl)) #cluster standard errors



#10 km
#Distance to Mysore-Bombay Border
rd5.mb=data[which(data$NEAR_DIST_border1<5000),] 

table(rd5.mb$border1)

#Distance to Hyderabad-Bombay Border
rd5.hb=data[which(data$NEAR_DIST_border2<5000),] 

table(rd5.hb$border2)


#Mysore
#outcome-high schools
educ.mys5=lm(high_binary~border1+TOT_POP+
              TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd5.mb) #OLS estimation
summary(educ.mys5)
educ.mys5.cl=cluster.vcov(educ.mys5, rd5.mb$dist_name)
educ.mys5.se=sqrt(diag(educ.mys5.cl)) #cluster standard errors

#Hyderabad
#outcome-high schools
educ.hyd5=lm(high_binary~border2+TOT_POP+
              TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd5.hb) #ols estimation 
summary(educ.hyd5)
educ.hyd5.cl=cluster.vcov(educ.hyd5, rd5.hb$dist_name)
educ.hyd5.se=sqrt(diag(educ.hyd5.cl)) #cluster standard errors




#50 km
#Distance to Mysore-Bombay Border
rd25.mb=data[which(data$NEAR_DIST_border1<25000),] 

table(rd25.mb$border1)

#Distance to Hyderabad-Bombay Border
rd25.hb=data[which(data$NEAR_DIST_border2<25000),] 

table(rd25.hb$border2)

#Mysore
#outcome-high schools
educ.mys25=lm(high_binary~border1+TOT_POP+
                TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd25.mb) #OLS estimation
summary(educ.mys25)
educ.mys25.cl=cluster.vcov(educ.mys25, rd25.mb$dist_name)
educ.mys25.se=sqrt(diag(educ.mys25.cl)) #cluster standard errors



#Hyderabad
#outcome-high schools
educ.hyd25=lm(high_binary~border2+TOT_POP+
                TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd25.hb) #ols estimation 
summary(educ.hyd25)
educ.hyd25.cl=cluster.vcov(educ.hyd25, rd25.hb$dist_name)
educ.hyd25.se=sqrt(diag(educ.hyd25.cl)) #cluster standard errors


#100 km


#Distance to Mysore-Bombay Border
rd50.mb=data[which(data$NEAR_DIST_border1<50000),] #20 km

table(rd50.mb$border1)

#Distance to Hyderabad-Bombay Border
rd50.hb=data[which(data$NEAR_DIST_border2<50000),] #20 km

table(rd50.hb$border2)


#baseline bandwidth
#Mysore
#outcome-high schools
educ.mys50=lm(high_binary~border1+TOT_POP+
                TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd50.mb) #OLS estimation
summary(educ.mys50)
educ.mys50.cl=cluster.vcov(educ.mys50, rd50.mb$dist_name)
educ.mys50.se=sqrt(diag(educ.mys50.cl)) #cluster standard errors



#Hyderabad
#outcome-high schools
educ.hyd50=lm(high_binary~border2+TOT_POP+
                TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd50.hb) #ols estimation 
summary(educ.hyd50)
educ.hyd50.cl=cluster.vcov(educ.hyd50, rd50.hb$dist_name)
educ.hyd50.se=sqrt(diag(educ.hyd50.cl)) #cluster standard errors







stargazer(educ.mys5, educ.mys10, educ.mys25, educ.mys50, 
          educ.hyd5, educ.hyd10, educ.hyd25, educ.hyd50,
          se=list(educ.mys5.se, educ.mys10.se, educ.mys25.se, educ.mys50.se,
            educ.hyd5.se, educ.hyd10.se, educ.hyd25.se, educ.hyd50.se), digits=3,
          omit=c("TOT_POP", "TOT_SC", "TOT_ST", "Slope", "TerrainRug", "Latitude", "Longitude"),
          dep.var.labels="Education: High Schools Availability",
          column.labels =c("bw:10km", "bw:20km", "bw:50km", "bw:100km",
                           "bw:10km", "bw:20km", "bw:50km", "bw:100km"),
          covariate.labels = c("Indirect Rule (Mysore)", "Indirect Rule (Hyderabad)", "Constant"),
          add.lines = list(c("Controls","\\checkmark","\\checkmark","\\checkmark","\\checkmark",
                             "\\checkmark","\\checkmark","\\checkmark","\\checkmark")),
          omit.stat = c("rsq", "f", "adj.rsq", "ser"))





